Models Cholesterol_prior_sens and
asprin_logistic_random_1_priors_sensitivity originally
failed with Attempt to redefine node errors. I got them to
run by moving the troublesome lines to data blocks. These chunks are
marked with eval=TRUE.
library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots') # caterplot
library('dplyr')
asprin_normal_fixed_1.txtThis is the model from Figure 4.1 (with a few extra digits of precision in the data).
model_code <- "model
{
for (i in 1:Nstud)
{
P[i] <- 1/V[i]
y[i] ~ dnorm(d, P[i])
}
d ~ dnorm(0, 1.0E-5)
OR <- exp(d)
}" %>% textConnection
# DATA
data <- list(
y=c(.3289011, .3845458, .2195622, .2222206, .2254672, -.1246363, .1109658),
V=c(.0388957,.0411673,.0204915,.0647646,.0351996,.0096167,.0015062),
Nstud=7
)
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
list(d=0)
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 7
## Unobserved stochastic nodes: 1
## Total graph size: 27
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50% 75%
## OR 1.1154942 0.03697400 1.04527705 1.09023746 1.1149566 1.1404389
## d 0.1087484 0.03314299 0.04428197 0.08639553 0.1088155 0.1314132
## deviance -3.7395092 1.41023676 -4.73853030 -4.63900994 -4.2773680 -3.4146738
## 97.5% Rhat n.eff overlap0 f
## OR 1.189913 NA 1 0 1.0000
## d 0.173880 NA 1 0 0.9996
## deviance 0.288420 NA 1 1 0.9704
plot(results)
asprin_normal_random_1.txtThis is the model from Figure 4.3, except it uses lower case for the y[] variable vector.
model_code <- "model
{
for (i in 1:Nstud)
{
P[i] <- 1/V[i]
y[i] ~ dnorm(delta[i], P[i])
delta[i] ~ dnorm(d, prec)
}
d ~ dnorm(0, 1.0E-5)
OR <- exp(d)
tau ~ dunif(0,10)
tau.sq <- tau*tau
prec <- 1/(tau.sq)
}" %>% textConnection
# DATA (the same as for the fixed effects model above)
data <- list(
y=c(.3289011, .3845458, .2195622, .2222206, .2254672, -.1246363, .1109658),
V=c(.0388957,.0411673,.0204915,.0647646,.0351996,.0096167,.0015062),
Nstud=7
)
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
list(d = 0, tau = 1, delta = c(0,0,0,0,0,0,0))
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d", "prec", "tau.sq"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 7
## Unobserved stochastic nodes: 9
## Total graph size: 38
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR 1.160938e+00 1.122747e-01 9.791676e-01 1.091816793 1.14621597
## d 1.447896e-01 9.327927e-02 -2.105249e-02 0.087843091 0.13646606
## prec 1.199015e+04 5.899705e+05 5.368759e+00 22.382149023 51.48911368
## tau.sq 3.712591e-02 5.701788e-02 7.957887e-05 0.006451449 0.01942158
## deviance -7.852934e+00 3.351394e+00 -1.306269e+01 -10.498470035 -8.27432876
## 75% 97.5% Rhat n.eff overlap0 f
## OR 1.21432738 1.4237089 NA 1 0 1.0000000
## d 0.19419032 0.3532654 NA 1 1 0.9621333
## prec 155.00394007 12566.1511858 NA 1 0 1.0000000
## tau.sq 0.04467846 0.1862628 NA 1 0 1.0000000
## deviance -5.48017083 -0.7749141 NA 1 0 0.9838000
plot(results)
asprin_logistic_random_1.txtThis is the model from Figure 4.5 (except that the line
OR <- exp(d) is in a different place).
model_code <- "model
{
for( i in 1 : Nstud ) {
rA[i] ~ dbin(pA[i], nA[i])
rB[i] ~ dbin(pB[i], nB[i])
logit(pA[i]) <- mu[i]
logit(pB[i]) <- mu[i] + delta[i]
mu[i] ~ dnorm(0.0,1.0E-5)
delta[i] ~ dnorm(d, prec)
}
OR <- exp(d) # !!! this line comes last in Figure 4.5
d ~ dnorm(0.0,1.0E-6)
tau ~ dunif(0,10)
tau.sq <- tau*tau
prec <- 1/(tau.sq)
}" %>% textConnection
# DATA
data <- list(
rA=c(49,44,102,32,85,246,1570),
rB=c(67,64,126,38,52,219,1720),
nA=c(615,758,832,317,810,2267,8587),
nB=c(624,771,850,309,406,2257,8600),
Nstud=7
)
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
list(d=0, tau=1, delta=c(0,0,0,0,0,0,0), mu=c(0,0,0,0,0,0,0))
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d", "prec", "tau.sq"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 16
## Total graph size: 74
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR 1.166672e+00 1.187119e-01 0.9813849004 1.095177e+00 1.15035462
## d 1.493962e-01 9.621368e-02 -0.0187905424 9.091595e-02 0.14007026
## prec 1.492055e+04 5.250138e+05 5.4773630417 2.195923e+01 48.64527803
## tau.sq 3.808501e-02 6.082000e-02 0.0001689156 7.508704e-03 0.02055698
## deviance 1.044143e+02 5.056343e+00 96.3293298236 1.007238e+02 103.84752807
## 75% 97.5% Rhat n.eff overlap0 f
## OR 1.22119445 1.4374136 NA 1 0 1.0000000
## d 0.19982943 0.3628454 NA 1 1 0.9635333
## prec 133.17880903 5920.1166164 NA 1 0 1.0000000
## tau.sq 0.04553893 0.1825696 NA 1 0 1.0000000
## deviance 107.44671303 115.9898590 NA 1 0 1.0000000
plot(results)
asprin_logistic_random1_predict.txtmodel_code <- "model
{
for( i in 1 : Nstud ) {
rA[i] ~ dbin(pA[i], nA[i])
rB[i] ~ dbin(pB[i], nB[i])
logit(pA[i]) <- mu[i]
logit(pB[i]) <- mu[i] + delta[i]
mu[i] ~ dnorm(0.0,1.0E-5)
delta[i] ~ dnorm(d, prec)
}
OR <- exp(d)
d ~ dnorm(0.0,1.0E-6)
tau~dunif(0,10)
tau.sq<-tau*tau
prec<-1/(tau.sq)
d.new ~ dnorm(d, prec)
OR.new <-exp(d.new)
}" %>% textConnection
# DATA
data <- list(
rA=c(49,44,102,32,85,246,1570),
rB=c(67,64,126,38,52,219,1720),
nA=c(615,758,832,317,810,2267,8587),
nB=c(624,771,850,309,406,2257,8600),
Nstud=7
)
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
list(d=0, tau =1, delta=c(0,0,0,0,0,0,0), mu=c(0,0,0,0,0,0,0))
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "OR.new", "d.new", "delta", "tau.sq"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 17
## Total graph size: 76
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR 1.16482673 0.11727540 0.978232230 1.093407115 1.15076623
## OR.new 1.18844654 0.30154239 0.737259168 1.040355891 1.14289962
## d.new 0.14666203 0.22323556 -0.304815798 0.039562857 0.13356856
## delta[1] 0.20910410 0.14090445 -0.031117108 0.111941872 0.19040952
## delta[2] 0.22763201 0.14769461 -0.015417530 0.122630883 0.20781802
## delta[3] 0.17825449 0.11126175 -0.019969587 0.102755952 0.16876713
## delta[4] 0.16310130 0.14905091 -0.119012662 0.072870934 0.15053163
## delta[5] 0.17095064 0.12942249 -0.070930538 0.087314807 0.15946323
## delta[6] -0.02882947 0.10186744 -0.238741923 -0.098195039 -0.02536030
## delta[7] 0.11311926 0.03778339 0.037838460 0.087995759 0.11337927
## tau.sq 0.04009716 0.06411076 0.000241616 0.008282165 0.02212963
## deviance 104.33057397 4.97948307 96.211258608 100.750212169 103.83358607
## 75% 97.5% Rhat n.eff overlap0 f
## OR 1.21896888 1.4333353 NA 1 0 1.0000000
## OR.new 1.27782385 1.8985090 NA 1 0 1.0000000
## d.new 0.24515851 0.6410688 NA 1 1 0.8200000
## delta[1] 0.29308486 0.5260593 NA 1 1 0.9574000
## delta[2] 0.31808150 0.5573990 NA 1 1 0.9666667
## delta[3] 0.24675171 0.4168797 NA 1 1 0.9593333
## delta[4] 0.24749593 0.4904861 NA 1 1 0.8915333
## delta[5] 0.24855173 0.4524951 NA 1 1 0.9276667
## delta[6] 0.04576045 0.1522363 NA 1 1 0.5945333
## delta[7] 0.13863053 0.1865140 NA 1 0 0.9986667
## tau.sq 0.04677193 0.1944284 NA 1 0 1.0000000
## deviance 107.32915712 115.4040223 NA 1 0 1.0000000
plot(results)
asprin_logistic_random_1_priors_sensitivityRUNTIME ERROR:
Compilation error on line 13.
Attempt to redefine node rBx[1,1]
I moved the offending lines into a data block, as recommended on stackoverflow.
Original WinBUGS code:
model {
## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
for (j in 1:3)
{
for( i in 1 : Nstud )
{
rAx[j,i] ~ dbin(pA[j,i], nA[i])
rAx[j,i] <- rA[i]
rBx[j,i] ~ dbin(pB[j,i], nB[i])
rBx[j,i] <- rB[i]
logit(pA[j,i]) <- mu[j,i]
logit(pB[j,i]) <- mu[j,i] + delta[j,i]
mu[j,i] ~ dnorm(0.0,1.0E-5)
delta[j,i] ~ dnorm(d[j], prec[j])
}
}
for (j in 1:3)
{
OR[j] <- exp(d[j])
d[j] ~ dnorm(0.0,1.0E-6)
m[j]~dnorm(0.0,1.0E-6)
}
# Setting three different priors for the between study variance
prec[1] <- 1/tau.sq[1]
tau.sq[1] <- tau[1]*tau[1]
tau[1] ~ dnorm(0,1.0E-6)I(0,)
prec[2] ~ dgamma(0.001,0.001)
tau.sq[2] <- 1/prec[2]
tau[2] <- sqrt(tau.sq[2])
tau[3] ~ dunif(0,10)
tau.sq[3] <- tau[3]*tau[3]
prec[3] <- 1/tau.sq[3]
}
model_code <- "
data {
for (j in 1:3) {
for( i in 1 : Nstud ) {
rAx[j,i] <- rA[i]
rBx[j,i] <- rB[i]
}
}
}
model{
## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
for (j in 1:3) {
for( i in 1 : Nstud ) {
rAx[j,i] ~ dbin(pA[j,i], nA[i])
# rAx[j,i] <- rA[i]
rBx[j,i] ~ dbin(pB[j,i], nB[i])
# rBx[j,i] <- rB[i]
logit(pA[j,i]) <- mu[j,i]
logit(pB[j,i]) <- mu[j,i] + delta[j,i]
mu[j,i] ~ dnorm(0.0,1.0E-5)
delta[j,i] ~ dnorm(d[j], prec[j])
}
}
for (j in 1:3){
OR[j] <- exp(d[j])
d[j] ~ dnorm(0.0,1.0E-6)
m[j]~dnorm(0.0,1.0E-6)
}
# Setting three different priors for the between study variance
prec[1] <- 1/tau.sq[1]
tau.sq[1] <- tau[1]*tau[1]
tau[1] ~ dnorm(0,1.0E-6)I(0,)
prec[2] ~ dgamma(0.001,0.001)
tau.sq[2] <- 1/prec[2]
tau[2] <- sqrt(tau.sq[2])
tau[3] ~ dunif(0,10)
tau.sq[3] <- tau[3]*tau[3]
prec[3] <- 1/tau.sq[3]
}
" %>% textConnection
# DATA
data <- list(
rA=c(49,44,102,32,85,246,1570),
rB=c(67,64,126,38,52,219,1720),
nA=c(615,758,832,317,810,2267,8587),
nB=c(624,771,850,309,406,2257,8600),
Nstud=7
)
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
# WinBUGS code has extra closing parentheses (???)
list(
prec=c(NA,0.1,NA),
tau=c(0.01, NA, 0.01),
mu = structure(
.Data=c(0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0),
.Dim=c(3,7)),
delta = structure(
.Data=c(0,0,0,0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0,0,0),
.Dim=c(3,7)),
d =c(0.1,0.1,0.1)
)
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000
)
##
## Processing function input.......
##
## Done.
##
## Compiling data graph
## Resolving undeclared variables
## Allocating nodes
## Initializing
## Reading data back into data table
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 42
## Unobserved stochastic nodes: 51
## Total graph size: 200
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR[1] 1.1662222 0.11739451 0.981212714 1.09437789 1.149769
## OR[2] 1.1526115 0.09525176 0.997237644 1.09143046 1.141165
## OR[3] 1.1650587 0.11863985 0.976503633 1.09352852 1.149952
## d[1] 0.1489909 0.09671732 -0.018966010 0.09018606 0.139561
## d[2] 0.1387560 0.08024680 -0.002766179 0.08748919 0.132050
## d[3] 0.1479067 0.09766691 -0.023776808 0.08940964 0.139720
## deviance 313.4448436 8.64787962 297.834584100 307.38539084 312.924119
## 75% 97.5% Rhat n.eff overlap0 f
## OR[1] 1.2216751 1.4400013 NA 1 0 1.0000000
## OR[2] 1.2004897 1.3714749 NA 1 0 1.0000000
## OR[3] 1.2194936 1.4332097 NA 1 0 1.0000000
## d[1] 0.2002229 0.3646440 NA 1 1 0.9625333
## d[2] 0.1827296 0.3158868 NA 1 1 0.9730000
## d[3] 0.1984357 0.3599164 NA 1 1 0.9605000
## deviance 319.0349227 331.7744037 NA 1 0 1.0000000
plot(results)
chol_model.txt, chol_data.txtError parsing model file:
syntax error on line 5 near "var"
I’m not sure why this was in the Supplemental Materials. It seems to
be exactly the same as Cholesterol_fixed.txt below (Model
4.1 in the Exercises), except it uses the name ‘var’ instead of ‘V’. It
turns out that ‘var’ is problematic in JAGS; maybe it is a reserved word
or something? I changed it to V (as in
Cholesterol_fixed.txt), which seems to work.
model_code <- "model
{
for (i in 1:k)
{
precision[i] <- 1/V[i]
logor[i] ~ dnorm(theta, precision[i])
}
theta ~ dnorm(0, 1.0E-5)
OR <- exp(theta)
}
" %>% textConnection
# DATA
data <- list(k=34)
data_df <- "logor V
0.4893853 0.1987757
-0.0683811 0.054427
-0.4876372 0.0730886
-0.1630902 0.1108992
1.469676 1.186739
0.0057637 0.0608936
-0.1496961 0.0799553
-1.94591 2.351288
-0.0470209 0.0194635
-0.3824268 0.0552552
-0.5508629 0.0615112
-0.0984401 0.0510244
0.2939126 0.1708592
-0.4883528 1.582492
-0.4354064 0.0030637
-0.2885568 0.0496847
0.4696847 0.0716345
0.2692525 0.0101539
-1.563976 0.8819983
-0.3517397 0.3733933
-0.0480808 0.0298812
-1.109251 2.687944
0.0129744 0.0469846
-0.4124823 0.0383304
-0.0233855 0.0200697
0.0803459 0.0082181
-0.2820902 0.0420423
-1.098612 2.707904
0.521464 2.696409
-0.2501722 0.477032
1.025401 0.3643915
-0.7528253 0.0676265
-1.580711 2.272708
0.5030903 0.1426682" %>% read.delim(text=., sep="\t")
for (colname in names(data_df)){
data[[colname]] <- data_df[[colname]]
}
results <- jags(
data = data,
# inits = initial_values,
parameters.to.save = c("OR", "theta"),
model.file = model_code,
n.chains = 1, # length(initial_values),
# n.adapt = 100,
n.iter = 50000,
n.burnin = 20000 #, n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 34
## Unobserved stochastic nodes: 1
## Total graph size: 108
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50% 75%
## OR 0.8450758 0.02730002 0.7929252 0.8263773 0.8446826 0.8633376
## theta -0.1688505 0.03229762 -0.2320263 -0.1907038 -0.1687943 -0.1469495
## deviance 83.6143430 1.40159071 82.6188310 82.7223690 83.0754879 83.9443157
## 97.5% Rhat n.eff overlap0 f
## OR 0.8996653 NA 1 0 1
## theta -0.1057325 NA 1 0 1
## deviance 87.5968764 NA 1 0 1
plot(results)
Cholesterol_fixed.txtThis is the file for Model 4.1 from the Exercises, referenced on p. 89.
Model 1: Generic fixed effect model applied to the cholesterol lowering dataset
model_code <- "model
{
for (i in 1:k)
{
P[i] <- 1/V[i]
logor[i] ~ dnorm(d, P[i])
}
d ~ dnorm(0, 1.0E-5)
OR <- exp(d)
}" %>% textConnection
# DATA: re-used from previous model. Here the name of the second column in the data table is `V` instead of `var`. I changed that in the previous model because the name `var` seemed to be causing problems.
# data <- list(k=34)
#
# data_df <- "logor V
# 0.4893853 0.1987757
# -0.0683811 0.054427
# -0.4876372 0.0730886
# -0.1630902 0.1108992
# 1.469676 1.186739
# 0.0057637 0.0608936
# -0.1496961 0.0799553
# -1.94591 2.351288
# -0.0470209 0.0194635
# -0.3824268 0.0552552
# -0.5508629 0.0615112
# -0.0984401 0.0510244
# 0.2939126 0.1708592
# -0.4883528 1.582492
# -0.4354064 0.0030637
# -0.2885568 0.0496847
# 0.4696847 0.0716345
# 0.2692525 0.0101539
# -1.563976 0.8819983
# -0.3517397 0.3733933
# -0.0480808 0.0298812
# -1.109251 2.687944
# 0.0129744 0.0469846
# -0.4124823 0.0383304
# -0.0233855 0.0200697
# 0.0803459 0.0082181
# -0.2820902 0.0420423
# -1.098612 2.707904
# 0.521464 2.696409
# -0.2501722 0.477032
# 1.025401 0.3643915
# -0.7528253 0.0676265
# -1.580711 2.272708
# 0.5030903 0.1426682" %>% read.delim(text=., sep="\t")
#
# for (colname in names(data_df)){
# data[[colname]] <- data_df[[colname]]
# }
# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
list(d=0)
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 34
## Unobserved stochastic nodes: 1
## Total graph size: 108
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50% 75%
## OR 0.8447898 0.02718736 0.7929180 0.8262318 0.8446778 0.8628069
## d -0.1691850 0.03217537 -0.2320354 -0.1908800 -0.1688000 -0.1475643
## deviance 83.6067409 1.39355862 82.6187542 82.7180516 83.0653805 83.9429405
## 97.5% Rhat n.eff overlap0 f
## OR 0.8981821 NA 1 0 1
## d -0.1073824 NA 1 0 1
## deviance 87.3827943 NA 1 0 1
plot(results)
Cholesterol_random.txtThis is the file for Model 4.2 from the Exercises, referenced on p. 90.
Model 2: Generic random effect model applied to the cholesterol lowering dataset
model_code <- "model
{
for (i in 1:k)
{
P[i] <- 1/V[i]
logor[i] ~ dnorm(delta[i], P[i])
delta[i] ~ dnorm(d, prec)
}
d ~ dnorm(0, 1.0E-5)
OR <- exp(d)
tau ~ dunif(0,10)
tau.sq <- tau*tau
prec <- 1/(tau.sq)
d.new ~ dnorm(d, prec)
}" %>% textConnection
# DATA: re-used from previous chunks.
initial_values <- list(
list(d=0,
delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
d.new=0)
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 34
## Unobserved stochastic nodes: 37
## Total graph size: 147
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50% 75%
## OR 0.8912361 0.06032464 0.7768036 0.8508168 0.8892777 0.92939818
## d -0.1174281 0.06755136 -0.2525678 -0.1615585 -0.1173457 -0.07321802
## deviance 26.8779939 6.35016609 15.4871660 22.4505696 26.5334845 30.85735636
## 97.5% Rhat n.eff overlap0 f
## OR 1.0157291 NA 1 0 1.0000000
## d 0.0156067 NA 1 1 0.9590667
## deviance 40.3147275 NA 1 0 1.0000000
Cholesterol_random_logistic.txtThis is the file for Model 4.3 in the Exercises, p90.
model_code <- "model
{
for( i in 1 : k ) {
rA[i] ~ dbin(pA[i], nA[i])
rB[i] ~ dbin(pB[i], nB[i])
logit(pA[i]) <- mu[i]
logit(pB[i]) <- mu[i] + delta[i]
mu[i] ~ dnorm(0.0,1.0E-5)
delta[i] ~ dnorm(d, prec)
}
OR <- exp(d)
d ~ dnorm(0.0,1.0E-6)
tau~dunif(0,10)
tau.sq<-tau*tau
prec<-1/(tau.sq)
}
" %>% textConnection
# DATA
data <- list(k=34)
data_df <- "nB nA rB rA
50 50 17 12
285 147 70 38
156 119 37 40
123 129 20 24
54 26 8 1
427 143 81 27
199 194 28 31
30 33 0 3
424 422 174 178
206 206 41 55
244 253 31 51
350 367 42 48
47 48 23 20
23 29 1 2
5552 2789 1025 723
1149 1129 37 48
221 237 39 28
5331 5296 236 181
88 30 2 3
71 72 5 7
1906 1900 68 71
94 94 0 1
2051 2030 44 43
279 276 61 82
1018 1015 111 113
4541 4516 269 248
421 417 49 62
48 49 0 1
94 52 1 0
79 78 4 5
6582 1663 33 3
204 202 28 51
30 60 0 4
311 317 19 12" %>% read.delim(text=., sep="\t")
for (colname in names(data_df)){
data[[colname]] <- data_df[[colname]]
}
initial_values = list(
list(d=0,
tau=1,
delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), # rep(0, data$k)
mu=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "theta"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000 #, n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 68
## Unobserved stochastic nodes: 70
## Total graph size: 317
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR 0.8931273 0.06185544 0.7784156 0.8515065 0.8905472
## deviance 385.6625551 10.82747134 365.8803706 378.0923377 385.1371885
## 75% 97.5% Rhat n.eff overlap0 f
## OR 0.9320761 1.022215 NA 1 0 1
## deviance 392.6321371 408.210312 NA 1 0 1
plot(results)
Cholesterol_random_logistic_2Also referenced on p90.
Model 3 (Cont) Question 5: Direct OR model applied to the cholesterol lowering dataset
model_code <- "
model
{
for( i in 1 : k ) {
rA[i] ~ dbin(pA[i], nA[i])
rB[i] ~ dbin(pB[i], nB[i])
logit(pA[i]) <- mu[i]
logit(pB[i]) <- mu[i] + delta[i]
mu[i] ~ dnorm(0.0,1.0E-5)
delta[i] ~ dnorm(d, prec)
}
OR <- exp(d)
d ~ dnorm(0.0,1.0E-6)
tau~dunif(0,10)
tau.sq<-tau*tau
prec<-1/(tau.sq)
# Storing theta.pooled in delta[37] to add to caterpillar plot
delta[37] <- d
# Replicate from a predictive distribution for a future study
delta[39] ~ dnorm(d, prec)
}
" %>% textConnection
# DATA
data <- list(k=34)
data_df <- "nB nA rB rA
50 50 17 12
285 147 70 38
156 119 37 40
123 129 20 24
54 26 8 1
427 143 81 27
199 194 28 31
30 33 0 3
424 422 174 178
206 206 41 55
244 253 31 51
350 367 42 48
47 48 23 20
23 29 1 2
5552 2789 1025 723
1149 1129 37 48
221 237 39 28
5331 5296 236 181
88 30 2 3
71 72 5 7
1906 1900 68 71
94 94 0 1
2051 2030 44 43
279 276 61 82
1018 1015 111 113
4541 4516 269 248
421 417 49 62
48 49 0 1
94 52 1 0
79 78 4 5
6582 1663 33 3
204 202 28 51
30 60 0 4
311 317 19 12" %>% read.delim(text=., sep="\t")
for (colname in names(data_df)){
data[[colname]] <- data_df[[colname]]
}
initial_values = list(
list(
d=0,
tau=2,
delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,NA,NA,NA,NA,1),
mu=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "theta"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000 #, n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 68
## Unobserved stochastic nodes: 71
## Total graph size: 318
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50% 75%
## OR 0.8944183 0.06199873 0.7784412 0.853385 0.8918417 0.933044
## deviance 385.5913474 10.71695257 366.1612898 378.094530 385.0791317 392.534913
## 97.5% Rhat n.eff overlap0 f
## OR 1.026008 NA 1 0 1
## deviance 408.177713 NA 1 0 1
plot(results)
Cholesterol_prior_sensThis file is referenced for Model 4.4, p91.
Model 4: Direct OR model with built in sensitivity analysis to prior
RUNTIME ERROR:
Compilation error on line 14.
Attempt to redefine node rBx[1,1]
Original WinBUGS code:
model
{
## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
for (j in 1:3)
{
for( i in 1 : k )
{
rAx[j,i] ~ dbin(pA[j,i], nA[i])
rAx[j,i] <- rA[i]
rBx[j,i] ~ dbin(pB[j,i], nB[i])
rBx[j,i] <- rB[i]
logit(pA[j,i]) <- mu[j,i]
logit(pB[j,i]) <- mu[j,i] + delta[j,i]
mu[j,i] ~ dnorm(0.0,1.0E-5)
delta[j,i] ~ dnorm(d[j], prec[j])
}
}
for (j in 1:3)
{
OR[j] <- exp(d[j])
d[j] ~ dnorm(0.0,1.0E-6)
m[j] ~ dnorm(0.0,1.0E-6)
}
# Setting three different priors for the between study variance
prec[1] <-1/tau.sq[1]
tau.sq[1] <-tau[1]*tau[1]
tau[1] ~ dnorm(0,1.0E-6)I(0,)
prec[2] ~ dgamma(0.001,0.001)
tau.sq[2] <- 1/prec[2]
tau[2] <- sqrt(tau.sq[2])
tau[3] ~ dunif(0,50)
tau.sq[3] <- tau[3]*tau[3]
prec[3] <- 1/tau.sq[3]
}
model_code <- "
data{
for (j in 1:3)
{
for( i in 1 : k )
{
rAx[j,i] <- rA[i]
rBx[j,i] <- rB[i]
}
}
}
model
{
## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
for (j in 1:3)
{
for( i in 1 : k )
{
rAx[j,i] ~ dbin(pA[j,i], nA[i])
# rAx[j,i] <- rA[i] # !!! Attempt to redefine node
rBx[j,i] ~ dbin(pB[j,i], nB[i])
# rBx[j,i] <- rB[i] # !!! Attempt to redefine node
logit(pA[j,i]) <- mu[j,i]
logit(pB[j,i]) <- mu[j,i] + delta[j,i]
mu[j,i] ~ dnorm(0.0,1.0E-5)
delta[j,i] ~ dnorm(d[j], prec[j])
}
}
for (j in 1:3)
{
OR[j] <- exp(d[j])
d[j] ~ dnorm(0.0,1.0E-6)
m[j] ~ dnorm(0.0,1.0E-6)
}
# Setting three different priors for the between study variance
prec[1] <-1/tau.sq[1]
tau.sq[1] <-tau[1]*tau[1]
tau[1] ~ dnorm(0,1.0E-6)I(0,)
prec[2] ~ dgamma(0.001,0.001)
tau.sq[2] <- 1/prec[2]
tau[2] <- sqrt(tau.sq[2])
tau[3] ~ dunif(0,50)
tau.sq[3] <- tau[3]*tau[3]
prec[3] <- 1/tau.sq[3]
}" %>% textConnection
# DATA
data <- list(k=34)
data_df <- "nB nA rB rA
50 50 17 12
285 147 70 38
156 119 37 40
123 129 20 24
54 26 8 1
427 143 81 27
199 194 28 31
30 33 0 3
424 422 174 178
206 206 41 55
244 253 31 51
350 367 42 48
47 48 23 20
23 29 1 2
5552 2789 1025 723
1149 1129 37 48
221 237 39 28
5331 5296 236 181
88 30 2 3
71 72 5 7
1906 1900 68 71
94 94 0 1
2051 2030 44 43
279 276 61 82
1018 1015 111 113
4541 4516 269 248
421 417 49 62
48 49 0 1
94 52 1 0
79 78 4 5
6582 1663 33 3
204 202 28 51
30 60 0 4
311 317 19 12" %>% read.delim(text=., sep="\t")
for (colname in names(data_df)){
data[[colname]] <- data_df[[colname]]
}
initial_values = list(
list(
prec=c(NA,0.1,NA),
tau=c(0.5, NA, 0.5),
mu = structure(.Data=c(0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0),.Dim=c(3,34)),
delta = structure(.Data=c(
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,
0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0, 0,0,0,0),.Dim=c(3,34)) #)
,
d =c(0.1,0.1,0.1)
)
)
results <- jags(
data = data,
inits = initial_values,
parameters.to.save = c("OR", "d", "theta"),
model.file = model_code,
n.chains = length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000 #, n.thin = 2
)
##
## Processing function input.......
##
## Done.
##
## Compiling data graph
## Resolving undeclared variables
## Allocating nodes
## Initializing
## Reading data back into data table
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 204
## Unobserved stochastic nodes: 213
## Total graph size: 875
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
##
## Calculating statistics.......
##
## Done.
summary(results)
## mean sd 2.5% 25% 50%
## OR[1] 0.8942644 0.06234855 0.7779298 0.8527763 0.8915863
## OR[2] 0.8927113 0.05906542 0.7816431 0.8530278 0.8910338
## OR[3] 0.8930435 0.06182829 0.7770039 0.8514910 0.8910776
## d[1] -0.1141740 0.06955541 -0.2511190 -0.1592580 -0.1147531
## d[2] -0.1156736 0.06604820 -0.2463570 -0.1589631 -0.1153729
## d[3] -0.1155074 0.06908539 -0.2523098 -0.1607664 -0.1153238
## deviance 1157.6300866 18.61739198 1122.9159571 1144.6421859 1157.0548490
## 75% 97.5% Rhat n.eff overlap0 f
## OR[1] 0.93307483 1.024227e+00 NA 1 0 1.0000000
## OR[2] 0.92989391 1.015263e+00 NA 1 0 1.0000000
## OR[3] 0.93218913 1.019841e+00 NA 1 0 1.0000000
## d[1] -0.06926988 2.393851e-02 NA 1 1 0.9502667
## d[2] -0.07268477 1.514736e-02 NA 1 1 0.9590667
## d[3] -0.07021956 1.964635e-02 NA 1 1 0.9543667
## deviance 1169.98851875 1.195988e+03 NA 1 0 1.0000000
plot(results)